# begin with the lists of PSB executives
executives<-read.table("061124-thesis-turnover.csv",header=T)

# Give all serving executives today's date
executives$E_D[executives$E_Y == "current"]<-format(Sys.time(),"%d")
executives$E_M[executives$E_Y == "current"]<-format(Sys.time(),"%m")
executives$E_Y[executives$E_Y == "current"]<-format(Sys.time(),"%Y")

# Assume that all executives for whom we lack month and day information
# were appointed on and removed on the 1st January
executives$S_D[is.na(executives$S_D)]<-1
executives$S_M[is.na(executives$S_M)]<-1

executives$E_D[is.na(executives$E_D)]<-1
executives$E_M[is.na(executives$E_M)]<-1

# join these together to create a date

executives$start<-as.Date(paste(executives$S_D,executives$S_M,executives$S_Y,sep="/"),format="%d/%m/%Y")
executives$end<-as.Date(paste(executives$E_D,executives$E_M,executives$E_Y,sep="/"),format="%d/%m/%Y")

# now move on to the government data set
# and repeat the same date operations

govts<-read.table("governments.csv",header=T,as.is=c(2:7))

govts$E_D[govts$E_Y == "current"]<-format(Sys.time(),"%d")
govts$E_M[govts$E_Y == "current"]<-format(Sys.time(),"%m")
govts$E_Y[govts$E_Y == "current"]<-format(Sys.time(),"%Y")

govts$S_D[is.na(govts$S_D)]<-1
govts$S_M[is.na(govts$S_M)]<-1

govts$E_D[is.na(govts$E_D)]<-1
govts$E_M[is.na(govts$E_M)]<-1

govts$start<-as.Date(paste(govts$S_D,govts$S_M,govts$S_Y,sep="/"),format="%d/%m/%Y")
govts$end<-as.Date(paste(govts$E_D,govts$E_M,govts$E_Y,sep="/"),format="%d/%m/%Y")

# begin to calculate indices
# first calculate TOR

# get the duration of each executive's term
# except for currently serving executives
executives$duration<-0
for(i in 1:nrow(executives)) {
	if(executives[i,"end"]==format(Sys.time(),"%Y-%m-%d")) (executives[i,"duration"]<-NA) else (executives[i,"duration"]<-executives[i,"end"]-executives[i,"start"])
}


# now get the average of these by country
# and divide it to get the rate of change every year

tor<-by(as.numeric(executives$duration),executives$country,mean,na.rm=T)

# now make a second pass to add long-serving current executives with tenures above average
# and recalculate tor

for(i in 1:nrow(executives)) {
	country<-executives[i,"country"]
	test1<-is.na(executives[i,"duration"])
	test2<-((executives[i,"end"]-executives[i,"start"])>tor[country])
	if(test1 & test2) (executives[i,"duration"]<-(executives[i,"end"]-executives[i,"start"]))
}

tor<-by(as.numeric(executives$duration),executives$country,mean,na.rm=T)
tor<-365/tor

# second calculate VUL
# create a notional boolean variable to see whether the government
# was shortly followed by a change in executive

govts$hitlist<-as.numeric(rep(0,nrow(govts)))

# now go through the list and see
# whether a change occured in the 180 days after govt format

for(i in 1:nrow(govts)) {
	test<-"NA"
	windowbegin<-govts[i,"start"]
	windowend<-windowbegin+180
	country<-as.character(govts[i,"Country"])
	test<-which(executives$start>windowbegin & executives$start<windowend & executives$country==country)
	if(length(test)==0) (govts[i,"hitlist"]<-0) else (govts[i,"hitlist"]<-1) 
}

vul<-by(govts[,"hitlist"],govts$Country,mean)

# you should now have two vectors with indice values

m.independence<-1-((tor+vul)/2)

# Read in other variables: Circ75, and PostCommunist

othervars<-read.table("othervars.csv",header=T)
# imputedvalue<-mean(othervars$PubSer,na.rm=T)
# othervars$PubSer[is.na(othervars$PubSer)]<-imputedvalue
# now look at DeJure

dejuretab<-read.table("dejure.csv",header=T)
dejuretab$aggopps<-rowMeans(dejuretab[,2:8],na.rm=T)
dejuretab$aggapps<-rowMeans(dejuretab[,9:14],na.rm=T)
# m.dejuremin<-apply(dejuretab[,15:16],MARGIN=1,FUN=min)
m.dejure<-rowMeans(dejuretab[,15:16])

# LaTeX output for the table of broadcasters and their independence scores
library(xtable)

table1<-cbind(as.vector(othervars$LongName),as.vector(othervars$Abbreviation),formatC(m.independence))
table1<-data.frame(table1)
table1<-table1[ order(m.independence,decreasing=T), ]
names(table1)<-c("PSB","Abbreviation","Independence")
print(xtable(table1,caption="PSBs listed by independence"),file="table1.tex",include.rownames=F,sanitize.text.function=function(x){x},caption.placement="top",size="scriptsize")

psbs<-cbind(othervars,as.vector(m.independence),m.dejure)
psbs<-data.frame(psbs)
names(psbs)<-c("country","Circ75","PostCommunist","HuIngPolariz93","PubSer","LongName","Abbreviation","Independence","DeJure")
attach(psbs)


# impute values
imputedvalues<-psbs
pubserreplacement<-mean(psbs$PubSer,na.rm=T)
imputedvalues$PubSer[is.na(imputedvalues$PubSer)]<-pubserreplacement
# now move on to the analysis
library(Zelig)
library(lmtest)
# full model
z.out1<-zelig(Independence~DeJure+log(Circ75)+log(Circ75):PostCommunist+HuIngPolariz93+PubSer,model="ls",data=psbs)
fullmodel<-summary(z.out1)
fullmodel
model1<-lm(psbs$Independence~psbs$DeJure+log(psbs$Circ75)+log(psbs$Circ75):psbs$PostCommunist+psbs$HuIngPolariz93+psbs$PubSer)
bptest(model1)

z.out2<-zelig(Independence~DeJure+log(Circ75)+log(Circ75):PostCommunist+HuIngPolariz93+PubSer,model="ls",data=imputedvalues)
fullmodel_imputed<-summary(z.out2)
fullmodel_imputed
model2<-lm(imputedvalues$Independence~imputedvalues$DeJure+log(imputedvalues$Circ75)+log(imputedvalues$Circ75):imputedvalues$PostCommunist+imputedvalues$HuIngPolariz93+imputedvalues$PubSer)
bptest(model2)

#reduced model
z.out3<-zelig(Independence~DeJure+log(Circ75)+log(Circ75):PostCommunist,model="ls",data=psbs)
model3<-lm(psbs$Independence~psbs$DeJure+log(psbs$Circ75)+log(psbs$Circ75):psbs$PostCommunist)
bptest(model3)

# Get the regression statistics
reducedmodel<-summary(z.out3)

# Start drawing table 2
# use parttable<-summary(z.out)
# parttable$coefficients[,c(1,4]

table2<-matrix(nc=6,nr=8)
table2[1:6,1:2]<-fullmodel$coefficients[,c(1,4)]
table2[1:3,4:5]<-reducedmodel$coefficients[1:3,c(1,4)]
table2[6,4:5]<-reducedmodel$coefficients[4,c(1,4)]
table2[8,1]<-fullmodel$adj.r.squared
table2[7,1]<-length(z.out1$residuals)
table2[8,4]<-reducedmodel$adj.r.squared
table2[7,4]<-length(z.out2$residuals)
table2<-data.frame(table2,row.names=c("Intercept","DeJure","logCirc75","Bureaucratic partisanship","Polarization","log(Circ75):PostCommunist","N","Adj. R-squared"))
names(table2)<-c("Coeff.","p. value","Sig.","Coeff.","p. value","Sig.")

table2[6,3]<-"**"
table2[2,6]<-"*"
table2[3,6]<-"**"
table2[6,6]<-"***"
xtable2<-xtable(table2,caption="Regression results",digits=3)
print(xtable2,file="table2.tex",include.rownames=T,include.colnames=T,caption.placement="top",hline=c(-1,-1,0,6,8,8))
# Start counterfactuals
# First counterfactual is for Gentiloni reform

x.out<-setx(z.out2,DeJure=0.66,PostCommunist=0,Circ75=142)
s.out<-sim(z.out2,x=x.out)
summary(s.out)

# Second counterfactual for Zapatero reforms

x.out<-setx(z.out2,DeJure=0.8,PostCommunist=0,Circ75=99)
s.out<-sim(z.out2,x=x.out)
summary(s.out)


# Appendix tables
tablea1<-data.frame(cbind(as.character(executives$Broadcaster),as.character(executives$Exec),as.character(executives$start),as.character(executives$end)))
names(tablea1)<-c("Broadcaster","Executive","Start","End")


xtablea1<-xtable(tablea1,caption="PSB executives")
print(xtablea1,file="tablea1.tex",include.rownames=F,include.colnames=T,caption.placement="top",floating=T,tabular.environment="longtable")
